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We revise the Levy’s construction of Brownian motion as a simple though still rigorous approach 
to operate with various Gaussian processes. A Brownian path is explicitly constructed as a linear 
combination of wavelet-based “geometrical features” at multiple length scales with random weights. 
Such a wavelet representation gives a closed formula mapping of the unit interval onto the functional 
space of Brownian paths. This formula elucidates many classical results about Brownian motion 
(e.g., non-differentiability of its path), providing intuitive feeling for non-mathematicians. The 
illustrative character of the wavelet representation, along with the simple structure of the underlying 
probability space, is different from the usual presentation of most classical textbooks. Similar 
concepts are discussed for fractional Brownian motion, Ornstein-Uhlenbeck process, Gaussian free 
field, and fractional Gaussian fields. Wavelet representations and dyadic decompositions form the 
basis of many highly efficient numerical methods to simulate Gaussian processes and fields, including 
Brownian motion and other diffusive processes in confining domains. 
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I. INTRODUCTION 

Diffusion is a fundamental transport mechanism in na¬ 
ture and industry, with applications ranging from physics 
to biology, chemistry, engineering, and economics. This 
process has attracted much attention during the last 
decades, particularly in statistical and condensed matter 
physics: diffusion-reaction processes; transport in porous 
media and biological tissues; trapping in heterogeneous 
systems; kinetic and aggregation phenomena like DLA, to 
name a few fields. From an intuitive point of view, Brow¬ 
nian motion is often considered as a continuous limit of 
lattice random walks. However, a more rigorous back¬ 
ground is needed to answer subtle questions. In math¬ 
ematical textbooks, Brownian motion is defined as an 
almost surely continuous process with independent nor¬ 
mally distributed increments [I|-|e|. The deceptive sim¬ 
plicity of this definition relies on the notion of “almost 
surely” that, in turn, requires a sophisticated formalism 
of Wiener measures in the space of continuous functions, 
filtrations, sigma-algebra, etc. Although this branch of 
mathematics is well developed, it is rather difficult for 
non-mathematicians, that is, the majority of scientists 
studying Brownian motion in their every-day research. 

In this paper, we discuss a different, but still rigor¬ 
ous, approach to define and operate with Brownian mo¬ 
tion as suggested by P. Levy (6j. We construct from 
scratch a simple and intuitively appealing representation 
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of this process that gives a closed formula mapping of 
the unit interval onto the functional space of Brownian 
paths. In this framework, sampling a Brownian path is 
nothing else than picking up uniformly a point from the 
unit interval. Figuratively speaking, Brownian motion is 
constructed here by adding randomly wavelet-based ge¬ 
ometrical features at multiple length scales. The explicit 
formula elucidates many classical results about Brown¬ 
ian motion (e.g., non-differentiability of its path). The 
illustrative character of the wavelet representation, along 
with the simple structure of the underlying probability 
space, is different from the usual presentation of most 
classical textbooks. 

Among various amazing properties, Brownian motion 
is known to have a self-similar structure: when a frag¬ 
ment of its path is magnified, it “looks” like the whole 
path. In other words, any fragment obeys the same prob¬ 
ability law as the whole path. As a consequence, Brown¬ 
ian paths exhibit their features at (infinitely) broad range 
of length scales. As a matter of fact, multiscale geo¬ 
metrical structures are ubiquitous in nature and mate¬ 
rial sciences [7j. For instance, respiratory and cardio¬ 
vascular systems start from large conduits (trachea and 
artery) that are then split into thinner and thinner chan¬ 
nels, up to the size of few hundred microns for the alve¬ 
oli and several microns for the smallest capillaries 
Another example is a high-performance concrete which 
is made with grains of different sizes (from centimeters 
to microns), smaller grains filling empty spaces between 
larger ones. The best adaptive description of such self¬ 
similar structures relies on intrinsicly multiscale func¬ 
tions to capture their mechanical or transport properties 
at different length scales. We illustrate this idea by con- 
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structing Brownian motion using wavelets, a family of 
functions with compact support and well defined scaling 
[MU. Wavelets appear as the natural mathematical lan¬ 
guage to describe and analyze multiscale structures, from 
heterogeneous rocks to biological tissues IT 2 MT 4 I I . The 
wavelet construction of Brownian motion naturally ex¬ 
tends to fractional Brownian motion and other Gaussian 
processes and fields, allowing one to efficiently simulate, 
for instance, turbulent diffusion with high Reynolds num¬ 
bers or financial markets. In particular, we discuss the 
Gaussian (or massless) free field and fractional Gaussian 
fields which appear as basic models in different areas 
of physics, from astrophysics (cosmic microwave back¬ 
ground) to critical phenomena, quantum physics, and 
turbulence [IM3- Written in a spectral form in one 
dimension, Brownian motion and the fractional Gaus¬ 
sian field look very similar, one of them being the frac¬ 
tional derivative of the other. Putting together these two 
processes reveals deep relations between them, and this 
correspondence carries over to higher dimensions. 

Most importantly, the wavelet representation is a start¬ 
ing point for a number of highly efficient numerical meth¬ 
ods to simulate various Gaussian processes and fields. 
Though wavelet representations and the related numeri¬ 
cal methods are all known, they are not easily available 
in a single source. Indeed the totality of these meth¬ 
ods seems to be poorly understood, even amongst spe¬ 
cialists. The purpose of this article is to present a uni¬ 
fied and intuitive framework that is based on elemen¬ 
tary mathematical structures like, for example, dyadic 
subdivision or wavelet tree. We also present some sim¬ 
ple number-theoretic shortcuts and consequent numerical 
algorithms. Finally, we discuss fast simulations of Brow¬ 
nian motion in confining domains, where one of the diffi¬ 
culties is the ability to quickly access the local geometry 
near the boundary. These techniques can be applied for 
studying Brownian motion and related processes or solv¬ 
ing partial differential equations in complex multiscale 
media. 

We hope that this didactic article will provoke inter¬ 
esting discussions amongst the experts and will help for 
a better understanding of both theory and implementa¬ 
tion of Brownian motion and other Gaussian processes 
and fields for a much broader community of their practi¬ 
cal “users”, namely, physicists, biologists, chemists, en¬ 
gineers, and economists. 


II. BROWNIAN MOTION 

In this section, we derive a wavelet representation of 
Brownian motion in a simple explicit way, allowing one 
to gain an intuitive feeling of this constructive approach. 


A. Physical view: upscaling and downscaling 

In order to illustrate the basic idea of a wavelet rep¬ 
resentation, we revisit the first single particle tracking 
experiment by R. Brown who looked through a micro¬ 
scope at stochastic trajectories (now known as Brownian 
paths) of pollens of Clarkia (primrose family) 0. First 
examples of such trajectories for mastic grains in water 
were reported by J. Perrin [20|, HI|. The essence of the 
wavelet representation can be recognized in his descrip¬ 
tion of three trajectories (Fig. [[]) that were recorded at 
30-second intervals [2lj : u Ils ne donnent qu’une idee tres 
affaiblie du prodigieux enchevetrement de la trajectoire 
reelle. Si, en effet, on faisait des pointes de seconde en 
seconde, chacun de ces segments rectilignes se trouverait 
remplace par un contour polygonal de 30 cotes relative- 
ment aussi complique que le dessin ici reproduit, et ainsi 
de suite." [ 96 | 

Following this idea, let us record the positions of a 
particle (e.g., pollen or grain) at successive time mo¬ 
ments with a selected time resolution 5. Each particle 
submerged in water is permanently “bombarded” by wa¬ 
ter molecules. Since the number of the surrounding water 
molecules is very large and their tiny actions are mostly 
uncorrelated, net microscopic displacements of the par¬ 
ticle cannot be considered deterministicly as in classical 
mechanics, but random. Since the interaction of water 
molecules between them is very rapid as compared to the 
macroscopic resolution scale, there is no memory effect 
in their action on the heavy particle. As a consequence, 
the microscopic displacements of the particle are (almost) 
independent and have a finite variance cr 2 . The macro¬ 
scopic displacement during the resolution time S is the 
sum of a large number N of these displacements with 
zero mean (no coherent flow). Although the average dis¬ 
placement is also zero, the stochastic fluctuations around 
this value are of the order of cr\fiN. Moreover, the central 
limit theorem gives us a precise probabilistic description 
of the fluctuations, resulting in a normal (or Gaussian) 
distribution of macroscopic displacements of the particle 
[22| . It is worth stressing that the Gaussian character of 
the macroscopic displacements appears without any spe¬ 
cific knowledge about the microscopic interactions. The 
only important information at microscopic level was sta¬ 
tionary, uncorrelated character of the interaction, and 
finite variance (if one of these conditions is missing, the 
resulting macroscopic process may exhibit anomalous be¬ 
havior, see f23|-[25} and references therein). This is known 
as coarsening or upscaling: complex interactions and the 
specific features of the underlying microscopic dynamics 
are averaged out on macroscopic scales. This is the rea¬ 
son why Brownian motion (or diffusion in general) is so 
ubiquitous in nature and science. Once we know that 
the microscopic details are irrelevant (under the condi¬ 
tions mentioned above), we can extend the Gaussian be¬ 
havior from macroscopic scales, where it has been es¬ 
tablished, to microscopic scales. This procedure can be 
called downscaling, when we explicitly and purposefully 
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FIG. 1: Three random trajectories of small mastic grains in 
water recorded by J. Perrin at 30-second intervals (reproduced 
from Hi]). 

transpose the universal macroscopic behavior even onto 
smaller scales. The resulting model of the microscopic dy¬ 
namics exhibits Gaussian features at all scales. While the 
true dynamics and its Gaussian model can be completely 
different at microscopic scales, they become identical at 
macroscopic scales. 

Knowing that a particle moves continuously, we con¬ 
nect its successive positions separated on time S by a 
continuous line. Since the experimental setup is limited 
to the selected time scale 5 , nothing can be said about 
the trajectory of the particle in between two records. In 
other words, the only condition for the trajectory to pass 
through the recorded points leaves us a variety of choices 
for the shape of the connecting continuous line. The com¬ 
mon choice is connecting the successive positions by lin¬ 
ear segments. As one will see, this choice fixes a partic¬ 
ular wavelet representation, the Haar wavelets. We shall 
show that other wavelets, corresponding to other choices 
of continuous connections, are as well useful. 

When the magnification and time resolution of the ex¬ 
perimental setup are increased, smaller details of the par¬ 
ticle’s trajectory appear, allowing one to refine the above 
piecewise linear approximation. Repeating this proce¬ 
dure, in theory up to infinity, one recovers all geometrical 
features and thus constructs the whole Brownian path. 
In what follows, we put this schematic description into a 
more rigorous mathematical frame. 


B. Mathematical view: multiscale construction 

We start with the position “records” at every unit 
time: t = 1,2,3,... (e.g., every second) along one co¬ 
ordinate (two- and higher dimensional Brownian motion 
is then obtained by taking d independent copies of the 
one-dimensional process). We focus on the time inter¬ 
val between 0 and 1, the construction being applicable 


to any interval [£, £+ 1]. For convenience, Brownian mo¬ 
tion is started at t = 0 from the origin: B( 0) = 0. By 
definition (or as a consequence of the central limit the¬ 
orem if one relies on the above physical reasoning), the 
position at time t = 1, S(l), is a random variable a o dis¬ 
tributed according to the standard normal (or Gaussian) 
law, AC(0,1), with mean zero and variance one 

e — z 2 /2 

P{ao £ ( x,x + dx )} = dx — _ _ . ( 1 ) 

v27T 

In physics, the variance a 2 is related to the diffusion co¬ 
efficient, D = a 2 /(2r), where r is the one step duration; 
here, cr = r = 1 yielding D = 1/2 throughout the pa¬ 
per. A linear approximation at this time scale (5 = 1) is 
simply 

B 0 (t) = a 0 t, 

that connects the positions B( 0) = 0 and -B(l) = ao by 
a linear segment. 

If the time resolution is doubled, a new, intermediate 
position b' = B( 1/2) can now be seen (Fig. [2]). The ran¬ 
dom variable b' is conditioned by the fact that Brownian 
motion passes through the points (0,0) and (l,ao), the 
value of ao being already known. It is distributed accord¬ 
ing to the normal law with mean value ^(B(0) + B(1)) = 
cio/2 and variance 1/4 (see Appendix lAl. In other words, 
one can write b' = ao/2 + aoo/2, where the new random 
variable aoo £ A/"(0,1) (i.e., distributed according to the 
standard normal law ([1|)) is independent of B( 0) =0 and 
B(l) = a 0 . 

The linear approximation at the time scale 6 = 1/2 
connects three successive points (0,0), (1/2,6') and 
(1, ao) by two linear segments: 

Mb', 0 <t<§, 

2t(ao — b') + (2b' — ao), \ < t < 1. 

The shape of this approximation looks like a skewed tent 
(Fig. [2]) that can be represented as a sum of a linear shift 
and a “symmetric tent” function: 

Bx(t) = a 0 t + (2b 1 - a 0 )h 00 (t) = a 0 t + a 00 h 0 o(t ), (3) 

where the “symmetric tent” function hoo(t) is 

it, 0<t<±, 

hoo(t) = O - 1, \<t< 1, (4) 

I 0, otherwise. 

The decomposition (J3J into the linear function t and the 
tent function hoo(t) is unique. The new approximation 
B\(t) is obtained from the previous one, Bo(t), by adding 
the new term representing a smaller geometrical detail. 

The same concept is applicable at every scale. As¬ 
sume that an approximation B n (t ) of Brownian motion 
is already constructed at the time scale S = 2~ n , i.e., 
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FIG. 2: Iterative construction of Brownian motion: (a) lin¬ 
ear approximation Bo(t) by a segment at scale 2°; (b) linear 
approximation Bi(t) by two segments at scale 2 _1 . The lat¬ 
ter “skewed tent” can be uniquely represented as the sum of 
a “symmetric tent” and a linear shift. 

the positions b k = B{t k ) are known at successive times 
ffc = k2~ n , k ranging from 0 to 2 n . The approximate 
Brownian path is a piecewise linear function passing suc¬ 
cessively through these points. 

At the next time scale 2 _n_1 , a new, intermediate 
position b' = B(t') of Brownian motion at time t' = 
{ ffc + ffc+i)/2 = (k + l/2)2~ n should be determined for 
each k. As previously, the random variable b' is condi¬ 
tioned by the fact that Brownian motion is known to pass 
through the points ( tk,b k ) and {t k +i, 6/c+i). It is again 
distributed according to the normal law, with mean value 
{bk + bk+ 1)/2 and variance 2 ™/4. In other words, one 
can write b' as 

b' = — {bk + bk+\) + 2 " I 2 1 a n k , (5) 

where the new normal random variable a n k € A/”(0,1) is 
independent of the other positions. The linear segment 
between {tk, bk) and (tfc+i, bk+i) is then replaced by two 
linear segments connecting the three successive points 
( tk,bk), ( t',b '), and (tk+i, b k +i)- This is a new “skewed 
tent” function which can be uniquely decomposed as the 
sum of the previous linear segment (drift) and a symmet¬ 
ric tent function h n k{t ) with the weight a n k , where 

h nk {t) = 2- n ^h 00 {2 n t-k) (6) 

is a rescaled symmetric tent function on the interval 
Ink = [k2~ n , {k + l)2 -n ) (see Fig. [dji). We stress again 
that the new approximation is obtained from the previous 
one by simply adding the tent function h nk (t), represent¬ 
ing a smaller geometrical detail at the new scale 2~ n ~ 1 , 


FIG. 3: Tent function h n k{t) (a) and the related Haar func¬ 
tion H n k{t) (b), both having the support on the interval 
Ink = [k2~ n , {k + 1)2 _ "). 


weighted by a normally distributed coefficient a n k which 
is independent of the previously determined positions. 

This construction is applicable to all linear segments 
(all k) at the given scale n , and it is valid for any scale. 
Repeating this procedure from the scale 2° up to infinity, 
one obtains the Haar wavelet representation of Brownian 
motion on the unit interval: 

oo 2" —1 

B“{t)=a“ 0 t + Y J Y.<M t )i ( ? ) 

n =0 k —0 

where all weights cLq and a“ k are independent Af{0, 1) 
random variables. Here we introduced the superscript ui 
in order to stress that a sampling of Gaussian weights a nk 
yields a random realization of Brownian motion. We will 
discuss in Sect. Ill Dl that all these Gaussian weights can 
be constructed from a single random number ui from the 
unit interval that provides a natural parameterization of 
Brownian paths. 

The dyadic structure of the intervals implies that for 
any n £ Z, there exists only one interval I nk of length 
2~ n containing the point t: k2~ n < t < {k + l)2~ n . 
The index k is simply the integer part of 2 n t: k = L2"tJ 
(i.e., the largest integer that does not exceed 2 n t). As 
a consequence, the convergence in the above formula is 
very rapid. In fact, if one needs to obtain the value of 
function B u (t) with a desired precision e, it is sufficient 
to calculate the first log 2 (l/e) terms, log 2 x being the 
logarithm of x on the base of 2. 

Subtracting the linear term a^t from Eq. © yields 
the Haar wavelet representation of a Brownian bridge on 
the unit interval, i.e., Brownian motion conditioned to 
return to 0 at time t = 1. 


C. Dyadic decomposition and interval subdivision 

In the wavelet representation ©, the first sum is car¬ 
ried over all scales n, while the second sum covers all 2 n 
subintervals I nk at the scale n. In some cases, it is con¬ 
venient to enumerate all the dyadic subintervals I nk by 
a single index m = 2 n + k as shown on Fig. 0] Since k is 
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FIG. 5: Generation of an infinite sequence of independent 
uniformly distributed variables ui 2 , u> 3 , a> 5 , • ■ • using binary 
expansion of a single number ui from the unit interval. For 
instance, 0 J 2 is constituted of the 2 nd , 4 th , 8 th , ... bits of oj. 


FIG. 4: Enumeration of the dyadic subintervals by a single 
index m. 


ranging from 0 to 2 ra_1 , the new index m uniquely iden¬ 
tifies the interval I n k . In particular, one easily retrieves 
the pair (n, k) from to as 

n = Ll°g2 m \ i k = to — 2™. 

Using the notations 

Og , to = 0 
®nfc. m >° 

we can write Eq. GD in a more compact form 

OO 

£“(t) = E < ~ h "*(*)■ ( 8 ) 

m =0 

As a result, Brownian trajectory is decomposed into a 
sum of tent functions (plus a linear drift) with random 
independent identically distributed Gaussian weights. As 
illustrated below, all these weights can be determined 
from a single uniformly distributed random variable. 



kmih) — 


t, TO = 0, 
hnkit), TO > 0, 


D. Representation of Gaussian weights 

In this subsection, we illustrate how all random Gaus¬ 
sian weights a“ can be explicitly related to a single uni¬ 
formly distributed random number 00 . In other words, 
we show that the complicated abstract probability space 
of Brownian paths can have a simple parameterization. 
However, this construction is not relevant from practical 
point of view, and the subsection can be skipped at first 
reading. 

We consider the binary expansion of a given real num¬ 
ber uj from the unit interval 


uniform picking up of oj in [ 0 , 1 ] is equivalent to inde¬ 
pendent random choice of its binary digits (or bits) bi. 
Then we choose a prime number p and construct another 
number uj p using the bits of oj at positions p, p 2 , p , ... 

uj p = Q.bpbp2 bp3 bpi ... ( 10 ) 

For example, ui 2 = 0 . 62 ^ 4 ^ 8^16 .. 0 J 3 = O .6369627681 • ■ 
etc. (Fig.[5D. If p and q are two different prime numbers 
then uip and ui q are independent as being constructed 
from separate sets of independent bits bi. Moreover, if oj 
is chosen uniformly from the unit interval, then each uj p 
is also uniformly distributed on the unit interval. Conse¬ 
quently, a single random number w gives rise to an infinite 
sequence {w p } of independent uniformly distributed ran¬ 
dom variables. In a more formal way, the real numbers 
< jj p can be written as 

OO 

cu p = ^ 2 -”- 1 (l + i?( 2 P' 1 u;)), ( 11 ) 

n= 1 

where R(x) = (—1)L X J is the Rademacher function. 

At last, we need to pass from uniformly distributed 
to normally distributed variables. For this purpose, we 
define the inverse $( 2 ;) of the error function: for x € 
[ 0 , 1 ], the value y of the function d>(a:) satisfies 


v 



— OO 


( 12 ) 


Although there is no simple analytic form for the func¬ 
tion $( 2 ;), many properties can be easily derived, and 
the whole function can be tabulated with any required 
precision. 

Ifp m denotes the (m+l) th prime (e.g., p 0 = 2 , pi = 3, 
P 2 = 5), then we set 


C = $ K™). ?n = 0,l,2,... (13) 


07 = 0 . 61626364 ... (9) 

where bi are equal to 0 or 1 (note that uj = 1 is expanded 
as 0.111... instead of its equivalent form 1.000 ...). A 


By construction, are independent normally dis¬ 
tributed random variables. In other words, Eqs. ED ED 
map a uniformly distributed oj onto a sequence of Gaus¬ 
sian weights a m . As a result, B UJ (t) is constructed as a 
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mapping from the unit interval, ui £ [0,1], onto the space 
of real-valued functions (more precisely, the Holder space 
H\/ 2 -e with any e > 0). In other words, any Brownian 
path is explicitly encoded by the real number ui. Picking 
up the real number to from the unit interval (with uniform 
measure) is thus equivalent to choosing a Brownian path 
(with Wiener measure). In this representation, the prob¬ 
ability space for Brownian motion is nothing else than the 
unit interval with uniform measure. It is intuitively much 
simpler than the classical construction of the probability 
space by means of Wiener measure, filtrations, etc. At 
the same time, this mapping is evidently neither continu¬ 
ous, nor injective (e.g., two numbers ui and u/ that differ 
at 6-th bit correspond to the same sequence of Gaussian 
random numbers). The mapping remains a rather for¬ 
mal construction whose only purpose was to show the 
equivalence between two spaces. 
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FIG. 6: Haar wavelets H n k(t) are obtained by dilations and 
translations of the mother function </> 1,1 (t) = Hoo(t) (shown 
on the left). 


E. Haar wavelets 


where dW u (t) denotes the Gaussian white noise which is 
defined here through the Haar wavelet decomposition 


In Eqs. © or ©, Brownian motion is decomposed into 
a sum of a linear function and tent functions. Figure [3]i 
illustrates that any tent function h nk {t) can actually be 
represented as the integral of a piecewise-constant func¬ 
tion 


t 

hnk{t) = J dt' H nk (t '), (14) 

o 

where H nk (t ) is called the Haar function and defined to 
be 0 on the complement of I n k and to take values 2 ra / 2 
and — 2"/ 2 on its left and right subintervals, respectively 
(Fig. ©:>). In fact, all Haar functions are obtained by 
translations and dilations of a single “mother” function 

[ 1, 0<t<i, 

= 1 , \<t< 1 , ( 15 ) 

i 0, otherwise. 

As illustrated on Fig. [6] one has 

H nk (t) = 2 n / 2 <f> 1 ' 1 (2 n t-k). (16) 

It is convenient to use previously introduced single in¬ 
dex to to denote different Haar functions: 


Hm(t) 


1, TO = 0, 
Hnk{t), TO > 0. 


Eq. © yields the following representation for Brownian 
motion 


t 



oo 

dW“(t) = t). (18) 

m —0 

It is easy to check that Haar functions {H m (t)} (to¬ 
gether with a constant function) form an orthonormal 
basis in the space L 2 ([0,1]) of measurable and square 
integrable functions that is 

l 

J dt (t) = . 

0 

Moreover, this basis is known to be complete in L 2 ([ 0,1]). 
This means that any function from L 2 ([ 0,1]) can be de¬ 
composed into a linear combination of the Haar functions 
(and a constant). From now on, we drop the superscript 
ui for getting simpler notations although the parameter¬ 
ization by ui remains valid for all discussed processes. 


F. General spectral representation 

The wavelet representation (fl8l) can be extended to 
any complete orthonormal basis {ipi(t)} of L 2 ([ 0,1]). In 
fact, if the orthonormal basis is complete, one can 

decompose any function H m (t) into a linear combination 
of 


OO 

Hm (t) = ^ ^ Cm,i 

i =0 

where the coefficients c mj i satisfy the orthogonality rela¬ 
tion 

oo 

= 1 (* = 0 , 1 , 2 ,...). 

m =0 


0 


(17) 


(19) 
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Substitution of this decomposition into Eq. m gives 

oo oo oo 

d,W(t) = ^2 dm^2c m ,i ipi(t) = (20) 

m =0 2=0 2 =0 

with new random weights 

OO 

di — ^ Cm,i- 

m— 0 

The sum of independent Gaussian variables is a Gaus¬ 
sian variable, and its variance is simply the sum of the 
squared coefficients i , which is equal to 1 according to 
Eq. (fl9l) . In other words, ai £ A/O,1). Moreover, the 
new random variables di are independent due to the or¬ 
thogonality of the functions We have thus shown 

that the Gaussian white noise can be decomposed into 
a linear combination with independent Gaussian weights 
in arbitrary complete orthonormal basis of L 2 ([ 0 , 1 ]). 

The completeness of the basis {4>i(t)} yields the usual 
covariance of the Gaussian white noise 

OO 

E{dW(ti)dW(t 2 )} = ^2 V’ii(*i)V’i 2 (* 2 )IE{aj 1 ai 2 } 
* 1 ,* 2=0 
OO 

= ^2 = S(ti - t 2 ), 

2 — 0 


where <5 is the Dirac distribution (or “^-function”). 
Substituting Eq. CHI) into Eq. 1171) . one gets 



dt' 4>i(t'), 


( 21 ) 


For instance, one can consider the Fourier basis on the 
unit interval, 


ipi{t) = a/2 cos(7r(z — 1/2)2), 

in order to get the Karhunen-Loeve expansion of Brow¬ 
nian motion .26]: 


£i < t 2 and £3 < £4 define two increments, B(t 2 ) — B(t\) 
and ^(£ 4 ) — B(t 3 ). If t 2 < £3 (i.e., the increments do not 
“overlay”), then they are independent (similar statement 
holds if £4 < £4 by symmetry). To proof this statement, 
we decompose the unit interval as 


[0,1] = [0, £ 1 ) U [£ 1 , £ 2 ] U (£ 2 ,£ 3 ) U [£ 3 , £ 4 ] U (£ 4 ,1], (23) 

(if one of subintervals [ 0 ,£ 1 ), (£ 2 ,£ 3 ) or (£ 4 , 1 ] is empty, 
it can be ignored). The basis {/>/£)} in Eq. (1211) can 
be chosen as the direct product of the Haar eigenbases 
on each subinterval. For instance, is the Haar 

basis of L 2 ([£i,£ 2 ]) which is extended to [0, 1]\[£i,£ 2 ] by 
zeros. In this particular representation, one has 


J. OO 

B(t 2 ) — B(t\) = / dt'^di 

ti i=0 

^ OO 

= / dt'j2^ i,t2] # >t2] (£') = 4 tl,t2] vt^Ti. 
1 2=0 


(24) 


In the second equality, the Haar functions from other 
subintervals (except [£i,£ 2 ]) vanished by construction. In 
turn, all the Haar functions ^f 1 ’* 2 ^') on [£i,£ 2 ] (with 
* > 0 ) vanished after integration due to their orthogonal¬ 
ity to the constant. The only remaining contribution is 
the constant term which has the unit L 2 ([£ 4 , £ 2 ]) norm: 
4^(t) = (t 2 — £i) -1 / 2 . Integrating this term, one gets 
the right-hand side of Eq. (l24l) which shows that the in¬ 
crement B(t 2 ) — B(t\) is a Gaussian variable with mean 
zero and variance £ 2 — £ 4 , as expected. Moreover, the 
same representation for [£ 3 ,£ 4 ] yields 


H(£ 4 ) - 5(£ 3 ) = i dt'^2 o-i V’i(i / ) 


*3 


2=0 


1 4 

/ OO 

d£'^di t3 ’‘ 41 ii2 3M] (t') = df M] 

is i=0 


(25) 


B(t) = 72 

2=0 


sin(7r(f — l/2)£) 
(i — 1 / 2 ) 7 t 


( 22 ) 


G. Basic properties of Brownian motion 


and the random weights Sq 1 ’* 2 ^ and Sq 3 ’* 4 ^ are indepen¬ 
dent by construction. As a consequence, the increments 
B(t 2 ) — B{t\) and B(ti) — B(t 3 ) are independent. 

(ii) The mean and covariance of Brownian motion are: 

E{B(£)} = 0, E{B(t!)B(t 2 )} = min {£ 4 , £ 2 }. (26) 


We have explicitly constructed the Haar wavelet de¬ 
composition 0 and then general representation m in 
order to reproduce the basic properties of Brownian mo¬ 
tion. Alternatively, one could first postulate such a rep¬ 
resentation as a definition of Brownian motion and then 
check that the basic properties are fulfilled. To illustrate 
this point, we check several properties. 

(i) Brownian motion is a Gaussian process with in¬ 
dependent increments. First, B(t) is Gaussian as be¬ 
ing a linear combination m of Gaussian variables. Let 


The first statement is obvious from Eq. (12T1) given that all 
weights have mean zero. To prove the second statement, 
one assumes that t 2 > £4 and considers 

E{5(£i)5(£ 2 )} = E {(Bih) - B(0))(B(t 2 ) - B(t 4 )} 

+ E{[H(£i)-f?(0)] 2 }. 

The first term vanishes due to the independence of incre¬ 
ments (and 5(0) = 0 ), while the second term is equal to 
£4 according to Eq. (El. 





(iii) Brownian motion is continuous but nowhere dif¬ 
ferentiable almost surely. The proof relies on a simple 
fact that the Gaussian weights a n k cannot be too large, 
e.g., the probability that |a n fc| > n decays extremely fast 
(as e~ n / 2 for large enough n). In turn, the norm of 
tent functions decreases exponentially that ensures the 
continuity of Brownian motion and the convergence of a 
partial sum approximation in Eq. 0 or similar expres¬ 
sions to Brownian motion. Moreover, the remainder of 
this approximation decreases exponentially fast with the 
truncated scale N (for technical details, see Appendix 

E). 


H. Alpert-Rokhlin wavelets 

As we mentioned in Sect. Ill Ai taking the particular 
orthonormal basis is equivalent to choosing a way to 
connect successive positions of Brownian motion at a fi¬ 
nite scale 6 . The Haar wavelets and the resulting tent 
functions present the simplest way of connection by lin¬ 
ear segments. Such a piecewise linear approximation of 
Brownian motion introduces singularities at the connec¬ 
tion nodes (corners). For some problems, it is convenient 
to deal with a smooth approximation of Brownian mo¬ 
tion at finite scales (although a true Brownian trajec¬ 
tory, the limiting curve, remains nowhere differentiable). 
For this purpose, one can use the Alpert-Rokhlin multi¬ 
wavelet basis (234301 - This basis is generated by a set 
of q functions </> 9,1 (t),..., cj> q,q (t ) which are supported on 
the interval [0,1], are piecewise polynomials of degree 
q— 1 on [0,1/2] and on [1/2,1], and satisfy the moment 
cancellation conditions 

1 

[ dt t k <fr q ' p (t) = 0, k = 0 ,l,...,q-l, ^ 

J w p=l,2,...,q. 

0 

These mother functions generate the Alpert-Rokhkin 
multiwavelets of order q by translations and dilations: 


lation conditions (l27l) : 




P’\t) 


■s/3 (1 — 4t), 0<f<±, 

\/3 (4i - 3), i < t < 1, 

0, otherwise 


6i — 1, 0 <t < \, 

6 t — 5, \ < t < 1, 

0, otherwise 


Higher-order mother functions (with q > 2) can be con¬ 
structed through an orthogonalization procedure (see 
[27l l3ll] for details and examples). 

Note that the wavelet representation of Brownian mo¬ 
tion involves the integral of wavelets 

t 

Kb P (t) = j dt’^ p (t'). (28) 

0 

For instance, one gets for q = 2 (Fig. 0 

v/3 i (1 — 2f), 0 <t<\, 

v / 3(l-t)(l-2t), \<t< 1, 

t{3t-l), 0<f<±, 

(t — l)(3t — 2), ±<t<l. 

while the other functions /i n ’ fc (t) and h% k (t) are obtained 
by dilations and translations. As a consequence, Brown¬ 
ian motion gets a closed formula in terms of the Alpert- 
Rohklin multiwavelets of order 2: 

00 2 n — 1 2 

B{t) = a 0 t + ai V3t(t-l) + Y^ J2 a nk h nk( t )’ ( 29 ) 

n =0 k— 0 p= 1 

where all weights do, a\, a^ k are independent A/"(0,1) 
variables, and the second term is the integral of the linear 
basis function \f3{2t — 1). An extension of this represen¬ 
tation to the Alpert-Rohklin multiwavelet basis of order 
q is straightforward. 


hoo (t) = 
h 2 oo(t) = 


cj ) q n P(t) = 2- n / 2 ^P(2 n t-k) 


n = 0,1,2,..., 
fc = 0,1,2,... 2" — 1. 


I. Numerical implementation 


The set of functions {^’^(t)}, completed by the set of 
orthonormal polynomials of order m < q, forms a com¬ 
plete basis of L 2 {[ 0,1]). This completion is necessary be¬ 
cause all mother functions (j) q,p {t ) (and thus all < Pnk(t )) 
are orthogonal by construction to all polynomials of or¬ 
der in < q. Similarly, the Haar wavelets were completed 
by a constant function. 

When q = 1, there is only one mother function </> 1,1 (t) 
defined by Eq. m which generates the Haar wavelets 
by translations and dilations. For q = 2, there are two 
mother functions (Fig. 0. satisfying the moment cancel- 


The wavelet representation of Brownian motion can be 
easily implemented in practice. To carry computations 
with a (fixed) desired precision e, it is sufficient to trun¬ 
cate the first sum in Eq. 0 or equivalent relation up 
to N = Llog 2 (l/e)J, because higher-order terms describe 
geometrical details at smaller scales. The remainder of 
this series can be estimated using Eq. (IB4I) . For Haar 
wavelets, such a truncation qualitatively corresponds to 
an approximation of Brownian motion by a broken line 
composed of linear segments of length close to e, while 
the Alpert-Rokhlin wavelets yield smoother approxima¬ 
tions (Fig. [0)- This is a fascinating feature of wavelets, 
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FIG. 7: Two mother functions for Alpert-Rokhlin first-order 
multiwavelets (on the left) and their integrals hgg’ ( t ) (on the 
right). The use of these wavelets corresponds to another type 
of connection (not by linear segments) between successive 
points of Brownian motion in the refinement procedure. For 
comparison, the tent function hoo(t) is shown on the last plot 
by dashed line. Note that horizontal and vertical scales are 
not matched here. 

allowing one to capture geometrical details at different 
scales. 

A realization of a Brownian path is completely de¬ 
termined by a set of random coefficients a“ fe (or a^J. 
In Sect. Eml we discussed an explicit scheme to gener¬ 
ate all these coefficients from a randomly chosen num¬ 
ber u> from the unit interval. In practice, one can use 
standard routines to generate pseudo-random normally 
distributed weights a“ fc (or ). The computation of 
tent functions h nk (t) can be easily implemented. Conse¬ 
quently, the computation of B(t) at any time t requires 
only log 2 (l/e) operations, each of them consisting of find¬ 
ing h nk (t), multiplying it by a“ fe , and summing their con¬ 
tributions. 

It is instructive to compare the wavelet approach to 
conventional techniques. We consider the computation of 
all positions b k = B{t k ) at equidistant times t k = k2~ N 
( k = 0,...,2 n ) at some scale e = 2~ N . In a classical 
scheme, Brownian motion is modeled by a sequence of 
small random jumps 

6 o = 0, b k+l =b k + 2 ~ N / 2 a' k (k = 0,1,..., 2^ — 1), 

( 3 °) 

with 2 n independent normally distributed random vari¬ 
ables a' k £ A/’(0,1). Similar computation relying on 
wavelet representations requires one random variable for 
a linear shift and 2 k random variables at each scale k, 
k ranging from 0 to IV — 1. The total number is then 
l + (l + 2-|-4-|-...-l-2 Ar_1 ) = 2 N . It is not surprising that 
both schemes require the same degree of randomness to 




FIG. 8: A random Brownian path at scales n = 3 (dashed 
line), and n = 10 (solid line) with Haar wavelets (a) and 
Alpert-Rohkin wavelets with q = 2 (b). 

represent a Brownian path at chosen scale. The wavelet 
representation does not reduce the complexity or ran¬ 
domness, but re-organize the data in a hierarchical struc¬ 
ture to facilitate their use. For instance, formula 0 ac¬ 
cesses approximate positions of Brownian motion at any 
time point t, not necessarily t k . In a classical scheme, one 
could use a linear interpolation between two neighboring 
points to get the same result. Again, the wavelets do not 
bring new features which are not available by conven¬ 
tional techniques, but provide another, structured and 
efficient, representation. 

Throughout the above sections, Brownian motion was 
constructed on the unit interval for convenience. The 
constructed process can be easily rescaled to any finite 
interval, while an extension to R+ or R. is possible as 
well. Finally, an extension to isotropic Brownian motion 
in R d is obtained by taking d independent samples of 
one-dimensional Brownian motion. 


III. BEYOND BROWNIAN MOTION 

A. Fractional Brownian motion 

A similar technique can be applied to construct and 
study fractional Brownian motion which is also known 
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as random fractal velocity field [H, [33|. For instance, 
random fractal velocity field with the Hurst exponent 
H = 1/3 (defined below) corresponds to the Kolmogorov 
spectrum in high Reynolds number turbulence [7], |3J-|36|. 
Fractional Brownian motion, as a Gaussian stochastic 
process with long-range correlations, has found numerous 
applications in different fields, ranging from transport 
phenomena in porous media 13 lL i37h 40| | to analysis of 
financial markets l4l| . 

P. Levy proposed the first extension of Eq. C3 by 
using the Riemann-Liouville fractional integration which 
can be thought of as a moving average of a Gaussian 
white noise 42] 


t 

B H {t) = r(g 1 + i J(t- s) H ~Uw(s), (31) 

where 0 < H < 1 is the Hurst exponent, and T(H + 1/2) 
is the normalization factor (r(z) being the Gamma func¬ 
tion). Mandelbrot and van Ness discussed the limitations 
of this definition (e.g., its strong emphasis on the origin) 
and proposed to use the Weyl fractional integral that 
yields j33| 


o 

Bh (*) = T ^ + i) | J ~ (-s) H ~^]dW{s) 

^ — OO 

t 

+ J(t- s) H ~?dW(s) 

0 

(32) 

for t > 0 (and similar for t < 0). The last representation 
can also be written as (see H3) 


t 

B H (t) = j K H (t, s) dW(s), (33) 

o 


where 


Knit , s) 


(t — s) H 2 

W+F 


2-Fl 


H— — , — 1 —- 

2 2 2 s 


(34) 


and 2 F 1 (a, 6; c; z) is the hypergeometric function. The 
ordinary Brownian motion is retrieved at H = 1/2 for all 
cases. 

Using a wavelet representation for the Gaussian 
white noise, one obtains 


B H (t) =2_^di I ds K H (t,s)4>i(s). 

*=0 n 


(35) 


The integrals can be evaluated using an appropriate ba¬ 
sis {tpi(t)}. Moreover, the moment cancellation property 
(1271) for the Alpert-Rokhlin multiwavelets with a large 


enough order q guarantees that the integrals in Eq. (l35l) 
are highly localized, yielding a rapid convergence of the 
above sum. This convergence is a key point for efficient 
numerical algorithms for simulation of fractional Brow¬ 
nian motion (see ;3l|, f37l - l39j l. Among other numerical 
methods, we mention alternative wavelet representations 
pl - l47j | (e.g., the method by Arby and Sellan is imple¬ 
mented in the Matlab function ‘wfrnb’), circulant em¬ 
bedding of the covariance matrix [48l - l50| . and random 
midpoint displacement method 5l| which is often used in 
computer graphics to generate random two-dimensional 
landscapes. 

In general, the kernel Kn{t^t') can be replaced by 
any convenient kernel to extend this approach to vari¬ 
ous Gaussian processes. For instance, setting K{t,t') = 
e -B(t-t ) yieidg ti ie Ornstein-Uhlenbeck process IaE [53{. 
Similarly, one can deal with various stochastic dynamics 
generated by Langevin equations [H]. 


B. Gaussian Free Field and its extensions 

Brownian motion is the integral of the Gaussian white 
noise which, in turn, is obtained as a linear combina¬ 
tion of orthonormal functions {tpi} forming a complete 
basis of the space L 2 ([0,1]), with standard Gaussian 
weights. This construction can be extended to any sepa¬ 
rable Hilbert space H. However, whatever the functional 
space H is taken, a linear combination of its orthonor¬ 
mal basis functions with standard Gaussian weights does 
not belong to this space (the argument is the same as for 
L 2 space, the norm of such a linear combination being 
infinite). In particular, the Gaussian white noise is not a 
function but a distribution. In order to construct exten¬ 
sions of the Gaussian white noise with desired properties, 
one needs to carefully choose the Hilbert space H. In this 
section, we briefly discuss two such extensions: Gaussian 
free field (GFF) {55] and fractional Gaussian field (FGF) 
with logarithmic correlations (5(| [57]]. The GFF ap¬ 
pears as the basic description of massless non-interacting 
particles in field theories. Both GFF and FGF consti¬ 
tute important models in different areas of physics, from 
astrophysics (describing stochastic anisotropy in cosmic 
microwave background) to critical phenomena, quantum 
physics, and turbulence [THI - ITsl . While the “sequence” of 
random variables of Brownian motion B (t) was naturally 
parameterized by “time” t (a real number from the unit 
interval or, in general, from R), random variables of a 
field can in general be parameterized by points from an 
Euclidean domain, a manifold, or a graph. For instance, 
one can speak about random surfaces which can model 
landscapes (e.g., mountains) or ocean’s water surface, in 
which the height is parameterized by two coordinates. 
The geometrical structure of “smooth” random surfaces 
was thoroughly investigated [58M60j | , especially for Gaus¬ 
sian fields which are fully characterized by their mean 
and covariance. Important examples of smooth Gaussian 
fields are the random plane waves and random spherical 
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harmonics (see [6l], fe). In turn, the GFF and FGF are 
examples of highly irregular random fields. As Brownian 
motion can be obtained as the limit of discrete random 
walks, the Gaussian free field in two dimensions appears 
in the limit of discrete random surfaces [gH HU • 

The Gaussian free field is constructed by choosing the 
Dirichlet Hilbert space , in which the scalar prod¬ 

uct of two functions / and g is defined as 

</,5>ffv(n) = J dx (V/-V.g) (36) 

a 

for a given Euclidean domain O C R d . When is 
bounded, an orthonormal basis of this space can be ob- 
tained by setting -0j( x) = 7 Ui(x ), where Ui{x) are the 

Abnormalized Dirichlet eigenfunctions of the Laplace op¬ 
erator: Aui(x) + XiUi(x) = 0 in Q with Ui(x) = 0 at the 
boundary dfl, and A i are the corresponding eigenvalues. 
The Gaussian free field on is defined as 

OO 

F(x) =^2 a i 4>i(x) (37) 

i—0 

where a.; £ A7(0,1) are independent Gaussian weights. 
Given that the eigenvalues A i asymptotically grow as i 2 / d 
according to the Weyl’s law [66l 1671 ], the sum in Eq. ED 
is convergent for d = 1 (in which case the GFF is simply 
a Brownian bridge) but diverges in higher dimensions. 
In the plane, this sum barely misses the convergence, be¬ 
ing logarithmically diverging. In quantum field theory, 
it is related to the infra-red divergence for massless par¬ 
ticles. As a consequence, F(x) is not a function but a 
distribution for d > 1. Being a linear combination of 
normal variables, the field F(x) is Gaussian and thus is 
fully characterized by its mean E{A(:r)} = 0 and the 
covariance 


E{F(xi)F(x 2 )} = yU A- 1 Ui(xi)ui(x 2 ) = G(xi,x 2 ), 
2=0 

( 38) 

where the right-hand side can be recognized as the Green 
function G(xi,x 2 ) of the Laplace operator in the domain 
O. Alternatively, one could define the GFF by setting the 
covariance equal to the Green function (in which case the 
definition holds even for unbounded domains). Strictly 
speaking, since the sum in Eq. ED diverges for d > 
1, the GFF should be treated as a distribution by its 
action on every fixed test function <j>{x). In particular, 
the covariance should be given by covariance of actions 
of F on two test functions 4>i and (f> 2 : 


(l36l) by fractional Laplacians 


</, g)HU Q) = Jdx ((-A)"/ • (-A Yg), (40) 

n 


for a positive v. The Dirichlet Hilbert space is re¬ 
trieved for v = 1/2. In a similar way, the functions 
tpi(x) = A ~ u ^ 2 Ui{x) form a basis of the space 
from which the FGF is constructed as in Eq. ED - This 
sum is convergent for v > d /4 and divergent otherwise. 
In the particular case v = d/ 4, the FGF is logarithmi¬ 
cally divergent in all dimensions. Note that the FGF 
coincides with the GFF in the plane (d = 2). In partic¬ 
ular, the FGF with logarithmic correlations on the unit 
interval reads explicitly as 


Fit) 



sin(7rfct) 

Vk 


(41) 


More generally, one can take any orthonormal ba¬ 
sis {iq(t)| of A 2 ([0,1]) and then apply the operator 

(—A) -1 / 4 = (^) to get the basis of AA^/ 4 ([0,1]). 
The advantage of the Laplacian eigenbasis is that the 
fractional integral operator (A) (which can be de¬ 
fined through the Fourier transform) is replaced by mul¬ 
tiplication by ^ 4 . As a consequence, the FGF with 
logarithmic corrections in one dimension appears as the 
half-derivative of Brownian motion: 



(42) 

This formula reveals a very close relation between these 
two processes which are often considered as distinct ob¬ 
jects. Note that the series in Eq. (TE1) diverges logarith¬ 
mically, while the derivative of order 1/2 — e would lead 
to converging series. In other words, Brownian motion 
belongs to the Holder space H 1 / 2 _ e (for any e > 0) so 
that its derivatives of order less than 1/2 exist, but 1/2 
and higher do not. As for Brownian motion, the explicit 
closed formula ED allows one to sample random realiza¬ 
tions of the FGF by picking up w from the unit interval, 
i.e., the probability space for this process is nothing else 
than the unit interval with uniform measure. 


IV. RESTRICTED DIFFUSION 


E{F(j)(xi), F(j)(x 2 )} = J dx\dx 2 G{xi—x 2 ) (j)i{xi) (j) 2 {x 2 ) 

f2xfi 

(39) 

The fractional Gaussian fields can be obtained by re¬ 
placing the gradient operators V in the scalar product 


In this section, we briefly discuss how multiscaling and 
dyadic decompositions help simulating restricted diffu¬ 
sion. This is an ubiquitous problem in physics (e.g., 
transport in porous media), chemistry (e.g., heteroge¬ 
neous catalysis), biology and physiology (e.g., diffusion 
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in cells, tissues, and organs). When Brownian motion is 
restricted, physico-chemical or biological interactions be¬ 
tween the diffusing particle and the interface of a confin¬ 
ing medium should be taken into account. For instance, 
paramagnetic impurities dispersed on a liquid/solid in¬ 
terface cause surface relaxation in nuclear magnetic res¬ 
onance (NMR) experiments jfisl . [iioj ; cellular membranes 
allow for a semi-permeable transport through the bound¬ 
ary (71)1472} ; chemical reaction may transform the particle 
or alter its diffusive properties [73 0- While an accurate 
description of these processes at the microscopic level is 
challenging, the contact with the interface is very rapid 
at macroscopic scales, allowing one to resort to an effec¬ 
tive description of the surface transport by an absorp¬ 
tion/reflection mechanism [zEEl- This mathematical 
process is known as (partially) reflected Brownian mo¬ 
tion IzzHzi- 

The presence of a boundary drastically changes the 
properties of Brownian motion (e.g., the reflecting 
boundary forces the process to remain inside the do¬ 
main) so that earlier wavelet representations cannot di¬ 
rectly incorporate the effect of the boundary. Since re¬ 
stricted diffusion is relevant for most physical, chemical 
and biological applications, various Monte Carlo meth¬ 
ods have been developed for simulating this stochastic 
process, computing the related statistics (e.g., the first 
passage times 0), and solving the underlying bound¬ 
ary value problems [8ll - l83 |. The slow convergence of 
Monte Carlo techniques (typically of the order of 1 /vM 
in the number of trials) requires fast generation of Brow¬ 
nian paths. The simplest generation of a Brownian path 
at successive times 6, 2 5, 3 5,... by adding normally dis¬ 
tributed displacements and checking the boundary effects 
at each step becomes inefficient in multiscale media. In 
fact, tiny geometrical details of the medium require the 
use of comparably small displacements, resulting in a 
very large number of steps needed to model large-scale 
excursions. 

To overcome this limitation, the concept of fast ran¬ 
dom walks was proposed [85]. The basic idea consists in 
adapting displacements to the local geometrical environ¬ 
ment, performing as large as possible displacements with¬ 
out violating the properties of Brownian motion. When 
the walker is at point x, the largest displacement is possi¬ 
ble at the distance \x — <9121 between x and the boundary 
<912 of an Euclidean domain O C R d . In fact, the ball 
B(x, \x — <9f2|) of radius \x — <9f2| does not contain any 
“obstacle” (e.g., piece of boundary) to the walker. Since 
Brownian motion is continuous, it must leave the ball be¬ 
fore approaching the boundary of the confining domain. 
The rotation symmetry implies that the exit points are 
distributed uniformly over the boundary of the ball. In¬ 
stead of modeling the fully-resolved trajectory of Brow¬ 
nian motion inside the ball, one can just pick up at ran¬ 
dom a point x' on the sphere of radius | | and move 

the random walker at this new position. The random 
duration of this displacement can be easily generated 
gaill. From here, one draws a new ball B(x ', \x' — <9f2|), 


and so on, until the walker approaches the boundary <912 
closer than a chosen threshold. From this point, an ap¬ 
propriate boundary effect (e.g., absorption, relaxation, 
chemical transformation, permeation, reflection, etc.) is 
implemented. Due to its efficiency, fast random walk 
algorithms have been used to simulate diffusion-limited 
aggregates (DLA) 86, 87f , to generate the harmonic mea¬ 
sure on fractals [73 (88} j89|, to model diffusion-reaction 
phenomena in spherical packs 0,m, to compute the 
signal attenuation in pulsed-gradient spin-echo experi¬ 
ments [92,[93j], etc. In this section, we focus on multiscale 
tools to estimate the distance, while other aspects of fast 
random walk algorithms can be found elsewhere [§4} . 


A. Distance to a boundary 

The efficiency of fast random walk algorithms fully re¬ 
lies on the ability to rapidly estimate the distance be¬ 
tween any point (e.g., the current position of the walker) 
and the boundary. Multiscale dyadic decompositions 
provide an efficient way to these estimates. To illustrate 
the idea, we first consider the one-dimensional case and 
then discuss its straightforward extension to the multidi¬ 
mensional case. 

In one dimension, the problem can be formulated as 
follows: given a set {x n } of N “boundary” points on the 
unit interval, how the distance to this set from another 
point x can be estimated in a rapid way? Successive 
computation of the distances \x — x n \ and finding their 
minimum is of course the simplest but the slowest way 
(of order of N). Instead of computing the distances to 
all boundary points, one can split the unit interval into 
two half-intervals, and check the distance to the points 
belonging to the half-interval that contains x. If the dis¬ 
tribution of points x n on [0,1] is more or less uniform, 
this division approximately halves the number of compu¬ 
tations. In the same spirit, splitting on subintervals of 
length 1/4, 1/8, etc. would reduce the number of compu¬ 
tations roughly by factors 4, 8, etc. Using such dyadic de¬ 
compositions, one needs approximately log 2 (iV) splitting 
to attend the level when one (or few) point x n belongs to 
the same subinterval as x. The number of computations 
is then of order of log 2 ( N ) (assuming the distribution of 
boundary points is more or less uniform). Moreover, if 
computation is carried with a desired precision e (to con¬ 
sider x n as “pointlike”, one needs s 1 /N ), one can con¬ 
tinue splitting up to the level log 2 (l/e) so that the length 
of subintervals becomes smaller than e. Since the points 
{x n } are stored with precision e, one cannot distinguish 
two points at any scale smaller than s. Consequently, 
any subinterval of length e can be either vacant, or oc¬ 
cupied by only one point x n (two points from the same 
subinterval would be indistinguishable). In this case, the 
number of computations, log 2 (l/e), is actually indepen¬ 
dent of whether the distribution of points x n is uniform 
or not. In other words, this algorithm can be applied for 
any finite set of points x n that are all distinguishable at 
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FIG. 9: Construction of a dyadic tree of subintervals for 
given boundary point x by its binary expansion. 




scale £, i.e., \x n — x m \ > £ for any n and m. 

B. Dyadic decomposition 

For practical implementation, the boundary points x n 
are used to generate a dyadic tree of subintervals at the 
scales ranging from 1 to log 2 (l/e). In turn, the binary 
expansion of the test point x is used to “navigate” search 
on the tree (Fig. 0). In fact, one can associate to a given 
point x G [0,1] a sequence of dyadic intervals I n ,[ 2 n xJ 
such that x £ I n ,\ 2 n x\ at any scale n. At each scale n, 
one chooses the left or the right subinterval depending on 
whether the n th bit is 0 or 1. Applying this procedure to 
all boundary points x n , one can generate a dyadic decom¬ 
position of the boundary. Figure II Oh shows an example 
with three boundary points {0, 0.4, 1} at five scales, 
from 2° to the smallest one e = 2~ 5 . The dyadic decom¬ 
position is stored as a tree, where a vertex is associated 
with a subinterval. Each vertex can be connected to one, 
two, or three other vertices (see Fig. IlOh ). The “height” 
of the vertex from the “root” determines the scale of the 
corresponding subinterval. 

Once a dyadic decomposition for a given boundary is 
constructed, it can be used to estimate the distance to the 
boundary from any point x. In fact, one can easily (and 
very rapidly) find the smallest subinterval I n \ 2 "x\ con¬ 
taining simultaneously x and some boundary point. For 
this purpose, one starts from the “root” vertex and de¬ 
scends on the tree using the bits of x to choose left or right 
edges at each scale. The descend is stopped when there is 
no edge to follow (Fig. [TUhb Once the smallest common 
interval I n ,\ 2 n x\ is found, there are two option: either 
n = |_log 2 (l/e)J so that the point x is indistinguishable 
from some boundary point at scale e, and the distance 
estimate is set to 0; or n < |_log 2 (l/e)J, and the distance 
to the boundary can be roughly estimated as 2~ n ~ 1 since 


FIG. 10: (a) Example of a dyadic decomposition of the 

unit interval with three boundary points {0, 0.4, 1} (shown 
by vertical dashed lines) and the related tree of subintervals 
at five levels (e = 2~ 5 ). For a test point x = 0.37 (shown 
by arrow), one uses its binary expansion x = 0 .01011... to 
navigate over the tree. The descend is stopped at level n = 2 
since x £ I33 = [3/8,4/8]. However, the rough estimate, 
2 - ™ -1 = 0.125, of the distance between x and the boundary, 
\x — <?fi| = 1 0.37 — 0.4| = 0.03, obviously fails. To apply the 
1/3-trick, one constructs the dyadic decomposition for the 
boundary points shifted by 1/3 (b). In this tree, the descent 
for the point x + 1/3 is stopped at level n' = 4. The combined 
lower estimate 2“ = 2“ 4 /6 « 0.0104 is valid. 


the next subdivision must separate the point x from the 
boundary points. However, this simplistic argument fails 
when the point x and the closest boundary point lie on 
opposite sides of the midpoint of the smallest common 
interval In,\ 2 n x\i but very close to each other. Although 
the next subdivision separates these two points, the dis¬ 
tance between them can be arbitrarily small so that the 
rough estimate 2~ n ~ 1 is wrong, as illustrated on Fig. 
m. To get the correct lower estimate, one can apply 
the so-called “1/3-trick”. 


C. The 1/3-trick 

The 1/3-trick can be easily illustrated for two points x 
and y. Let I n ,i 2 n x\ be the smallest common interval con¬ 
taining both points x and y. Suppose that the distance 
between these points is smaller than 2 -rl /6. We consider 
the points x' = x+ 1/3 and y' = y + 1/3 (shifted by 1/3), 
and determine their smallest common interval I n , y 2 n ' x'\ • 
As shown in AppendixjU] the distance between the points 
x' and y' (and thus between the points x and y) is larger 
than 2~ n /6. It is thus sufficient to find the smallest com- 
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mon intervals for the pair x, y and its shifted counterpart 
x',y', and the distance between the points is bounded 
below as 

\x — y\ = \x'-y'\ > 2- max 'b l .™'}/6. 

This simple fact allows one to rapidly estimate the dis¬ 
tance to boundary points. Let us consider a new bound¬ 
ary 0ft' which is obtained by shifting the old one by 1/3: 
dVL' = {ieR : i- 1/3 € 5ft}. For the new boundary, 
another dyadic tree of subintervals can be constructed 
(Fig. flOb l in order to estimate the distance between dCl' 
and the shifted point x' = x + 1/3 at a given scale e. 
While this construction may appear redundant at first 
thought because \x — dfl | = \x' — <9S7'|, the crucial point 
is that we search for a lower estimate at a finite scale 
at which two dyadic trees are different. Performing the 
descend over both dyadic trees (using the binary expa- 
sion of x and x ', respectively), one identifies the level n 
(resp. n') of the smallest common interval of x (resp. 
x') and the closest boundary point (resp. shifted closest 
boundary point). The lower estimate of the distance is 
then 

\x - 0ft| = \x' - 0ft'| > 2- max { n ’ n '>/6. 

D. Higher dimensions 

Similar constructions are applicable in higher dimen¬ 
sions which are more relevant for applications. The lower 
estimate relies on the generalized mean inequality: 

M q (c \,..., c n ) < M p (ci ,..., c n ) ( q<p ), (43) 

where the generalized mean M p (c \,..., c n ) of n positive 
numbers c\,... ,c n is 

/ i " \ 1/p 

M p (ci,...,c n ) = f ' ( 44 ) 

Setting q = 1, p = 2, and Ck = \xk — yu\ for two points 
x={xi ,...,x d ) and y = (yi,..., y d ), Eq. d43]) yields the 
lower estimate of their Euclidean distance: 

( d \ 1/2 i d 

\x-y\= ( - Uk) 2 I > -J= ^2 \x k - yk \• (45) 

Constructing two dyadic trees for each coordinate as 
earlier, one can then estimate each term \xk — yk\ and 
thus the distance from any point x = (xi,... ,x d ) to 
the boundary. In high dimensions (d 1), the right- 
hand side of the inequality (l45l) is strongly attenuated by 
the prefactor 1 /vd. This is an example of the so-called 
“curse of dimensionality”. To overcome this difficulty, 
one can implement random rotations and translations of 
the boundary. Although these transformations preserve 
the distance, they can improve the lower bound at a finite 


scale. Note that other multiscale constructions and the 
related searchable data structures can also be used such 
as Whitney decompositions of the computational domain 
(in which the size of each square (or cube) paving the 
domain is comparable to the distance to the boundary), 
quadtrees (or Q-trees), k-d trees, etc. {94j. 

E. Overall efficiency 

The advantages of multiscale dyadic trees are numer¬ 
ous: the simplicity of construction, the generality of 
boundary shapes, the rapidity of distance estimation, 
the flexibility for shape modifications, and low mem¬ 
ory usage. In fact, for a given precision e, the storage 
of the dyadic tree requires at worst N log 2 (l/e) inter¬ 
vals (i.e., log 2 (l/if) levels for each boundary point, for N 
points). In practice, this number is much smaller since 
many boundary points share the same interval at larger 
scales (e.g., the interval of size 1 is shared by all boundary 
points). 

The crucial point is that geometrical structure of the 
boundary does not matter at all: the method works for a 
random Cantor dust as well as for a circle. Moreover, the 
tree-like representation is highly adaptive, allowing one 
to modify the boundary from one set of simulations to 
other (or even from one run to the other). This feature 
can be very useful to study diffusion-controlled growth 
processes like DLA or transport phenomena in domains 
with moving boundaries. 

V. CONCLUSIONS 

In this paper, we revised the multiscale construction 
of Gaussian processes and fields. First, the Haar wavelet 
representation of Brownian motion was explicitly con¬ 
structed as a natural way to refine the geometrical fea¬ 
tures of a Brownian path under magnification. Since 
the Haar functions form a complete basis in the space 
L 2 ([0,1]) and their weights are Gaussian, such a repre¬ 
sentation can be extended to any complete basis {^»(t)} 
of L 2 ([ 0,1]), wavelet-like or not. In other words, the con¬ 
struction of Brownian motion has two separate “ingredi¬ 
ents”: deterministic functions capturing geometri¬ 

cal details, and their random weights ch. The choice of 
the functions ipi(t) is voluntary that gives a certain free¬ 
dom and flexibility in dealing with different problems. 
Qualitatively, this choice determines the way of connect¬ 
ing successive positions of Brownian motion at a given 
scale. On the opposite, the random weights determine 
the intrinsic stochastic properties of Brownian motion, 
independently of our choice of the basis 

The multiscale construction gives not only a simple 
closed formula for Brownian motion, but reveals its fun¬ 
damental properties. For instance, continuity and non¬ 
differentiability of Brownian paths naturally follow from 
this construction. In addition, we discussed a closed map- 
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ping from the sampling unit interval onto the space of 
Brownian paths. Sampling Brownian paths can there¬ 
fore be formally reduced to picking up a real number ui 
with the uniform measure. This construction does not re¬ 
quire elaborate notions from modern probability theory 
such as Wiener measures, sigma-algebras, filtrations, etc. 
Although these notions are useful, the explicit multiscale 
construction is much easier for non-mathematicians. 

These concepts are not limited to Brownian motion. 
We illustrated how fractional Brownian motion, Gaus¬ 
sian free field, and fractional Gaussian fields can be con¬ 
structed in a very similar way. Extensions to other Gaus¬ 
sian processes were also mentioned. Finally, we briefly 
discussed how the multiscale concepts can be used for 
simulating restricted diffusion (i.e., Brownian motion in 
confining domains). The dyadic subdivision and the re¬ 
lated hierarchical (multiscale) tree of subintervals allow 
one to rapidly estimate the distance to the boundary and 
thus to generate random displacements adapted to the 
local geometrical environment. Such fast random walk 
algorithms have found numerous applications. As for 
usual Brownian motion, dyadic decompositions appear 
as natural tools to store and rapidly access the geomet¬ 
rical information, resulting in fast algorithms. 
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Appendix A: Conditional law 

We explain why the distribution of the position of 
Brownian motion y = B( 1/2) at time t = 1/2 under con¬ 
ditions 5(0) = 0 and 5( 1) = x is given by the normal 
law Af(x/2, 1/4). 

Since the increments of Brownian motion on the unit 
intervals (0,1/2) and (1/2,1) are independent, the joint 
probability density p{B( 1/2) = y, B( 1) = x} is sim¬ 
ply equal to the product of the probability density 
p{B(1/2) — B(0) = y} to have the first increment equal to 
y and the probability density p{B(l) — B(l/2) = x—y} to 
have the second increment equal to x — y. These densities 
are given by normal laws with variance 1/2, yielding 


bility density we are looking for: 

(y-x/2) 2 /2 

P{B( 1/2) = y | 5(1) =x}= e E- ■ 

V27T Vl/4 

This is the normal distribution with the mean x/2 and 
the variance 1/4. 

Appendix B: Continuity and non-differentiability of 
Brownian motion 

In this section, we illustrate how the wavelet represen¬ 
tation of Brownian motion can be used to prove its basic 
properties such as continuity and non-differentiability. 

1. Continuity and convergence 

First, we show that a partial sum approximation con¬ 
verges to Brownian motion in both L°° and L 2 norms. 
Let us denote 

2 n — 1 

Fn(t) = £ (Bl) 

k =0 

so that B(t) = a 0 t + F 0 (t) + F 1 (t) + ... according to Eq. 
CD- Each function F n is a sum of hat functions with 
disjoint supports. To estimate the size of F n we need 
some upper bound on a n k■ We are going to show that 
\a n k\ < n for all sufficiently large n. 

Since a n k £ J\f(0, 1) are standard normal variables, we 
observe that 

— z 2 /2 

P(|a nfc | >n) = 2 J dz e ^_ < e~ n /2 

n 

for sufficiently large n. This inequality yields 

oo 2 n —1 oo 

y j y P(|a n fc| > n) < y 2 n e~ n2/2 < oo. 

n=0 k— 0 7i=0 

By Borel-Cantelli lemma this implies that |a n fc| < n for 
all but finitely many coefficients a n k- In other words, 
almost surely, there is finite, but random, N* such that 
\a n k\ < n if n > N*. For such n we have 

H^nlU-qo,!]) = 2 -n / 2-1 max{|a n fc|} < n2~ n / 2 ~ 1 . 

k 


p{B(l/2) = y, B(l)=x} 



( B — x / 2) 2 /2 

e ca 


V2^^/T74 


(x—y) 2 /2 \ 

e j 

V^y/172 J 

( e~* 2 / 2 \ 

)' 


Since the sum of these norms converges, the series ao t + 
^2 n F n (t) converges to B(t) in L°° norm (uniformly), i.e. 
for any e > 0 

^limPfllBti) - 5j V (t)ll L - ([0il]) > e} = °, (B2) 


N 2 n — 1 

B N (t) = a 0 t + EE ^nfchnfe(t) (B3) 

n=0 fe=0 


The second factor is simply the probability density for 
x = 5(1), while the first factor is the conditional proba- 


where 
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is a partial sum approximation of Brownian motion at 
scale 2~ N . This proves that Brownian motion is almost 
surely continuous. Moreover, the remainder of a partial 
sum approximation at scale 2 ~ N is exponentially small: 

\m)-B N m L ~< E IKWIL- < 

n=N +1 

(B4) 

when N is large enough. 

The L 2 convergence can be shown in the same way. 
We have 

2 n —1 

II^IIl 2 ([0,i]) = E l anfc | 2 H^ rafc llL 2 ([o,i]) < 2” n 2 2 “ n , 

fc =0 

where we used |a n fc| < n for large enough n. This in¬ 
equality implies the convergence of J7 II-P?i,||l 2 ([o,i]) with 
probability one and hence the series in Eq. 0 converges 
almost surely in L 2 norm: 


for all h. For a given scale n, let k be such that to is be¬ 
tween dyadic points t n ^~i and t n ^, where t n & = 2 ~ n k. 
The triangle inequality implies that for any j 

\B(t n ,k+j) ~ B(t n ,k+j-i)\ < \B{t n ,k+j) — B(t 0 )\ 

+ | B(t n , k+j -i) - B(t 0 )\ < M( 2 j + l) 2 ~ n . 

Let E n k be the event that this inequality holds for 
j = 1,2,3. Since increments are independent normal 
variables, one gets IP {E n k} < c( 2 -ra / 2 ) 3 , where c is a 
constant. The probability that these inequalities hold for 
some k from 0 to 2 n — 1 is then bounded by 2 n (c 2 ~ 3 n I 2 ) = 
c2 ~ n / 2 . The sum of these probabilities over n is finite, 
hence by Borel-Cantelli lemma, with probability one only 
finitely many of them will occur. On the other hand the 
assumed inequality (1B71) implies that infinitely many of 
E n k will occur. This yields the contradiction and proves 
that (EzD cannot be true. 

Appendix C: The 1/3-trick to estimate the distance 


jv lim o E{||B(t)- J B J v(t)|| i2([01]) } =0. (B5) 

Both statements (IB21IB5I) can be extended to arbitrary 
spectral representation of Brownian motion. Note also 
that these statements are applicable pointwise, e.g., 

lim E {\B(t) — -B/v(i)l} = 0 (B6) 

N-¥o O 


for any t € [0,1], 


2. Nowhere differentiability 


Since Brownian motion is a sum of hat functions it 
is easy to believe that B(t) it not differentiable almost 
everywhere. One can prove a much stronger statement 
that B{t) is nowhere differentiable with probability one. 
The proof goes along the same lines as for convergence 
(the argument follows the proof from (HI). 

We are going to show that almost surely for all t at 
least one of the two limits, 


b' (t) = lim sup 
B^(t) = lim inf 


B(t + h) - B(t ) 
h 

B{t + h) - B{t) 
h 


is infinite. This obviously implies that B{t) is almost 
surely nowhere differentiable. Note that at local extrema, 
one of these limits can be finite, so it is not true that both 
of them are always infinite. 

Let us assume that there is to such that both limits 
are finite at to ■ This implies that there is a random finite 
constant M such that 


B(to + h) — B{to) 
h 


(B7) 


Although the 1/3-trick is classical in analysis, we pro¬ 
vide some explanations which may be instructive for non¬ 
experts. 

The trick is based on a very simple result. Let I n k = 
[k2~ n , (k + l)2~ n ) and its boundary dl n k = {k2~ n , (k + 
l)2 -n }. For any real x and any integer n, there exist 
two intervals I n k and I n k> of the same length 2 ~ n that x 
belongs to I n k and x’ = x + 1/3 belongs to I n k’- If the 
distance from x to the boundary of I n k is smaller than 
2 -n /6, then the distance from x' to the boundary of I n k r 
is larger than 2 -n /6, and vice-versa. In other words, the 
point x and the shifted point x’ cannot be simultaneously 
close to the interval endpoints. 

Suppose the opposite is true, so that 

\2 n x-k\ < 1/6, |2"(x + 1/3) — k'\ < 1/6, 

where the integers k = k + [2(2 n x — k)] and k’ = k’ + 
[2(2 n x' — k')] denote the closest endpoint to x and x’, 
respectively. Since 2 n /3 = ko + (—l) n /3 with an integer 
ko , the point 2 n x should be simultaneously within the 
distance 1/6 to k and k' + ko + {— 1)"/3 that is impossible 
since the distance between these points is larger than 1/3 
(Fig. HU): 

i < \{k’ + k 0 + (-l) n /3) - k\ < \2 n x - k\ 

O 

+ \2 n x - {%’ + k 0 + (—1)73)| < 1/3- 

Using this simple result, one can prove the estimate 
for the distance from a given point x to the set of bound¬ 
ary points {y m }- Let I n k be the largest interval con¬ 
taining x and not containing any boundary point y m . 
Similarly, for the shifted point x' = x + 1/3, let I n ’k' be 
the largest interval containing x' and not containing any 
shifted boundary point y' m = y m + 1/3. Then the dis¬ 
tance from x to the set of boundary points {y m } is larger 
than 2 _max f"’ 7l, l/6. 


< M 
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2 n x 



\ 

k + 1/3 


The shifted boundary point y' belongs to some interval 
I n j. According to the previous result, the second inequal¬ 
ity implies that the shifted point y' cannot be close to the 
endpoints of I n j: 

2 _ "(i + 1/6) <y'< 2~ n (j + 5/6). 


FIG. 11: Illustration for the proof of the 1/3-trick. In this 
example, k = k' + ko = k and n is even. 

Suppose that n > n'. Assume that the statement is 
false so there exists a boundary point y € {y m } such 
that \x — y\ < 2“"/6. Suppose that y < x (the opposite 
case is similar). Since x € I n k and y </ I n k, both points 
x and y should be close to the endpoint k 2 ~ n : 

\x-k 2 - n \< 2 ~ n / 6 , \y-k 2 ~ n \< 2 - n / 6 . 


Since \x' — y'\ = |x — y\ < 2 n /6 and y < x, then y' < 
x' < y' + 2~ n /6, hence 

2 ~ n (j + 1/6) < x' < 2~ n (j + 1), 

i.e., the point x' belongs to the same interval I n j as y'. At 
the same time, x' belongs to I n >k' which is larger than 
I n j due to n > n!. The dyadic structure implies that 
Inj C I n 'k’ so that y' should belong to I n 'k' as well. But 
this is in contradiction with the initial assumption about 
In' k' • 
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